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HIGH-ORDER  FINITE  ELEMENT  METHODS  FOR 
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Abstract.  We  develop  a  framework  for  applying  high-order  finite  element  methods 
to  singularly-perturbed  elliptic  and  parabolic  differential  systems  that  utilizes  special  qua¬ 
drature  rules  to  confine  spurious  effects,  such  as  excess  diffusion  and  non-physical  oscilla¬ 
tions,  to  boundary  and  interior  layers.  This  approach  is  more  suited  for  use  with  adaptive 
mesh-refinement  and  order-variation  techniques  than  other  problem-dependent  methods. 
Quadrature  rules,  developed  for  two-point  convection-diffusion  and  reaction-diffusion 
problems,  are  used  with  finite  element  software  to  solve  examples  involving  ordinary  and 
partial  differential  equations.  Numerical  artifacts  are  confined  to  layers  for  all  combina¬ 
tions  of  meshes,  orders,  and  singular  perturbation  parameters  that  were  tested.  Radau  or 
Lobatto  quadrature  used  with  the  finite  element  method  to  solve,  respectively,  convection- 
and  reaction-diffusion  problems  provide  the  benefits  of  the  specialized  quadrature  formulas 

and  are  simpler  to  implement. 

Key  words,  adaptive  finite  element  methods,  Radau  and  Lobatto  quadrature,  high- 
order  methods,  singular  perturbations 
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1.  Introduction.  Beginning  with  a  trial  solution  of  a  differential  system  computed 
on  a  coarse  mesh  with  a  low-order  method,  an  adaptive  strategy  seeks  to  obtain  a  final 
solution  satisfying  prescribed  discretization  error  criteria  as  quickly  as  possible.  Tools  for 
“enriching”  the  initial  solution  are  denoted  as  h-,  p-,  or  r-refinement  when,  respectively, 
the  mesh  is  refined  and  coarsened,  the  order  of  the  method  is  varied,  or  the  mesh  is  moved 
to  follow  evolving  phenomena.  Enrichment,  typically  local  in  space,  can  be  either  local  or 
global  in  time  to  produce,  respectively,  a  local  refinement  method  [1]  or  a  method  of  lines 
[2].  The  most  successful  adaptive  enrichment  techniques  utilize  a  combination  of  the  base 
strategies  with  the  particular  combination  of  hp-refinement  capable  of  achieving  exponen¬ 
tial  convergence  rates  [2-7]. 

Singularly-perturbed  problems  are  ideal  differential  systems  for  adaptive  analysis 
because  it  is  far  more  efficient  to  resolve  the  nonuniform  behavior  within,  e.g.,  boundary 
layers  with  nonuniform  meshes  and  methods  than  with  uniform  structures.  Current  tools 
for  adaptive  analysis,  however,  lack  a  requisite  robustness  when  applied  to  singularly- 
perturbed  problems.  Symmetric  finite  difference  or  finite  element  methods,  commonly 
used  with  adaptive  software,  produce  spurious  oscillations  (§2)  when  applied  on  a  coarse 
mesh.  Based  on  this  incorrect  solution,  a  correct  error  estimation  could  indicate  global 
mesh  refinement.  While  possibly  successful,  this  strategy  is  far  from  optimal  since 
refinement  would  most  likely  be  necessary  only  within  layers.  “Upwinding”  eliminates 
oscillations  by  adding  diffusion,  but  successful  strategies  [8]  are  largely  unknown  for  other 
than  convection-diffusion  systems.  A  posteriori  error  estimates,  used  to  indicate  adaptive 
enrichment,  have  been  based  on  diffusive  dominance  for  both  elliptic  and  parabolic  sys¬ 
tems  [9-11];  thus,  minimally,  their  range  of  applicability  is  greatly  diminished  when  they 
are  applied  to  singularly-perturbed  problems. 

With  a  goal  of  developing  high-order  methods  that  may  be  used  with  adaptive  hp- 
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refinement  to  solve  singularly-perturbed  partial  differential  systems,  we  develop  a  Galerkin 
finite  element  technique  that  utilizes  special  quadrature  rules  to  attain  stability  (§2-4). 
Utilizing  a  hierarchical  framework  [12],  the  quadrature  rules  are  designed  to  integrate  pro¬ 
ducts  of  exponential  and  polynomial  functions  to  high  order.  In  this  regard,  they  bear 
some  relationship  to  exponentially-weighted  Petrov-Galerkin  methods  [13],  however,  poly¬ 
nomials  are  used  for  both  the  finite  element  trial  and  test  spaces.  Thus,  the  new  methods 
can  be  incorporated  within  an  existing  adaptive  finite  element  system  with  only  minor 
alterations.  The  singularly-perturbed  limit  of  the  derived  quadrature  rules  for  convection- 
diffusion  systems  are  Radau  integration  formulas  (§3).  Radau  (§3)  and  Lobatto  (§4)  qua¬ 
drature  yield  stable  computational  results  for,  respectively,  convection-  and  reaction- 
diffusion  problems  for  all  combinations  of  the  singularly-perturbed  parameter,  the  mesh 
spacing,  and  the  order  of  the  method.  Large  errors  are  confined  to  elements  containing  or 
adjacent  to  layers  and  may  be  reduced  by  either  h-  or  p-refinement.  Convergence  rates 
remain  optimal  in  the  diffusion  limit  and,  although  rates  are  unknown,  accuracy  is  high  in 
the  singularly-perturbed  limit. 


Quadrature  rules  are  derived  for  ordinary  differential  systems  and  our  methods  may 
be  useful  with  two-point  boundary  value  problem  software  [14].  However,  our  intent  is  to 
use  the  derived  rules  with  elliptic  and  parabolic  partial  differential  systems  through,  e.g., 
tensor  products  and  method  of  lines  reduction.  Thus,  after  verifying  that  the  specialized, 
Radau,  and  Lobatto  quadrature  rules  give  good  results  when  used  with  finite  element  pro¬ 
cedures  to  solve  ordinary  convection-diffusion  (§3)  and  reaction-diffusion  (§4)  systems, 
the  methods  are  applied  to  one-dimensional  transient  (§3)  and  two-dimensional  steady  (§5) 


problems.  Once  again,  the  results  display  a  robust  behavior  with  large  errors  confined  to 


regions  of  nonuniform  behavior. 

2.  Formulation  and  Background.  We  illustrate  the  quadrature-based  finite  element 


method  for  the  singularly-perturbed  two-point  boundary  value  problem 
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L u  :=  -eu"  +  c(x)u  +  d(x)u  =  f  (x),  a  <x  <b,  (2.1a) 

u(a)  =  uL,  u  (b)  =  uR ,  (2.1b, c) 

where  ( )'  :=  d{  )/dx.  Assume  that  e  >  0  and  consider  cases  when  (i)  c(x)  *  0  or  (ii) 
c(x)  :=  0  and  d(x)  >  0,  x  e  [ a,b ].  In  the  former  case,  convection  dominates  diffusion 
when  the  Reynolds  or  Peclet  number  ( b  -  a)  max  lc(x)/el  is  large,  and  the  solution  of 

(2.1)  has  a  boundary  layer  with  thickness  inversely  proportional  to  the  Peclet  number  at 
x  =  a  or  b  when  c  is,  respectively,  negative  or  positive.  Reaction  dominates  diffusion  in 
the  latter  case  when  ( b  —  a)[  max  ^(x)/£]^  is  large,  and  the  solution  of  (2.1)  features 

a<x<b 

boundary  layers  near  both  a  and  b . 

Consider  the  Galerkin  form  of  (2.1):  determine  u  e  Hg  satisfying 

B{y  ,u)  :=  £(vV)  +  ( v,cu  )  +  ( v,du )  =  ( vj ),  for  all  v  e  ,  (2.2a) 

where 

b 

(v,u)  =  Jv(x)u(x)dx,  (2.2b) 

a 

Hl  is  the  usual  Sobolev  space,  and  subscripts  £  or  0  restrict  functions  to,  respectively, 
satisfy  (2.1b,c)  or  trivial  versions  of  (2.1b,c).  A  finite  element  problem  is  constructed 
from  (2.2)  by  introducing  a  partition 

UN  :=  {  a  =  x0  <  xj  <•••<%=  b  }  (2.3) 

of  [a,  b]  into  N  subintervals;  defining  a  finite-dimensional  space 

SN,] p  =  { w  e  H 1  1  w(x)e  Tp.,  x  e  (x,-^,),  i  =  1,  2,  •••,  N  },  (2.4) 

where  p  :=  [px,  p2-  •••,  Pn  lr  and  Tp  is  a  space  of  polynomials  of  degree  p\  and  deter¬ 
mining  U  e  Se’p  satisfying 

B(V,U)  =  (V  J),  for  all  V  e  5^’p.  (2.5) 

Representing  SN  p  in  terms  of  a  hierarchical  basis  [12]  is  convenient,  efficient,  and 


stable.  Thus,  we  let 


where 


U(x)  =  £  Ci<bi(x)  cfoftx), 

i  =0  i—l  A:  =2 


$i(^(x)),  if  x  e  [*,■_!,*,•) 

4>/(x)  =  *  ^(^(x)),  if  ^  e  [JCi^i+i).  *  =  °* 

0,  otherwise 


(2.6a) 


(2.6b) 


J  if x  6  [*«-l*^i)  .  ..  c  \ 

<t>*(*)  =  |o,  otherwise  1  k=2,3,-,Pi,  *  =  >  ’  ’ 


&0  = 


2x  -  ATj _ !  -  X, 

Xi  -  X,_! 


$_i(S)  =  (i  - 1)/2,  $i(4)  =  a  +  4)/2, 


(2.7a) 


(2.7b, c) 


4©  =  ^  U,i^  -  5.  f-ui.  <™> 

and  Pfc(^)  is  the  Legendre  polynomial  of  degree  k.  The  piecewise  linear  basis  element 
<j>/(;c)  is  associated  with  the  node  xf  while  the  higher-degree  basis  elements  <t>*(*), 
k  _  2,  3,  •  ••,  pi,  are  associated  with  the  subinterval  (x,-^,).  In  order  to  simplify  the  sub¬ 
sequent  presentation,  let 

:=  {(J >/(x),  i  =0,  1  ,—,N,^(x),k  =  2,  3,  — , pt,  i  =  1,2 ,—,N}.  (2.8a) 

With  trivial  Dirichlet  boundary  conditions,  the  dimension  of  Sj?'p  is 


AT  =  -1  +  5>f. 

i= l 


(2.8b) 


Example  2.1.  Some  difficulties  that  arise  when  solving  (2.1)  by  a  conventional 
Galerkin-finite  element  approach  can  be  illustrated  for  the  simple  example  when  c  is  a 
constant,  d(x)  =  /(x)  :=  0,  UN  is  uniform  with  spacing  h,  and  pt  =  l,  i  =  1,  2,  N. 
In  this  case,  the  Galerkin  coordinates  satisfy 
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C0=UL,  CN=UR , 


(2.9a, b) 


cl+1  -  2c,  +  C,'.!  =  ■^'(CI+1  -  ci- 1)>  *  -  2’  ^ 


(2.9c) 


where  c,  =  t/ (*,),  i  =  0,  1,  N,  and 


ch 

P  =  T 


(2.9d) 


(2.9e,f) 


is  the  cell  Peclet  number.  The  solution  of  this  difference  equation  is 

V  -  XN  .  1  -  V  _  A  ,  ^  _  2  +  P 

7^'*’  ‘-0,1'  •*’  x‘2-p- 

The  computationally  challenging  singularly-perturbed  limit  occurs  when  p  »  1.  In  this 
case,  c,,  i  =  0,  1,  N,  oscillates  (i)  between  uL  and  uR  when  N  is  odd  or  (ii)  between 
the  line  joining  uL  and  uR  and  2p IN  when  N  is  even.  These  spurious  oscillations  are 
only  eliminated  when  the  mesh  is  sufficiently  fine,  i.e.,  when  p  is  reduced  to  0(1)  within 
the  boundary  layer. 


The  spurious  oscillations  can  be  eliminated  from  a  convection-diffusion  system  by 
introducing  a  directional  bias  that  adds  artificial  dissipation  to  the  system.  Although 
schemes  for  accomplishing  this  abound  for  convection-diffusion  f  systems  (cf.,  e.g.,  IT  in 
[15]),  Hemker’s  [13]  Petrov-Galerkin  formulation  has  a  generality  that  best  suits  our 
present  aims.  To  begin,  recall  that  the  Green’s  function  G(tjc)  for  the  operator  L  of  (2.1) 
is  an  element  of  H q  for  fixed  t  that  satisfies 

L *G  =  -eG^  -  [c(x)G]x  +  d[x)G  =0,  a  <x  <t,t  <x  <b,  (2.10a) 

[Gx(t,t)]x  =  f  :=  lim  [Gx(t,t+ 6)  -  Gx(t,t- 5)]  =  -j,  (2.10b) 

6— M3  £ 

where  an  x  subscript  denotes  partial  differentiation.  With  this  definition,  we  readily  show 
[13]  that 


v(r)  =  B(G(t,-),v),  for  all  veH1. 


(2.11) 
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Consider  the  Petrov-Galerkin  problem:  determine  U  e  5^  ,p  satisfying 

B(V,U)  =  (V,f),  for  all  V  e  Sq  ,p.  (2-12) 

Like  Sg’P,  the  space  So’p  is  a  finite-dimensional  subspace  of  Hq  ;  however,  it  consists  of 
piecewise  polynomial  and  exponential  functions.  We’ll  identify  a  basis  for  Sq  ,p  as 

*P  :=  {  vj//(jc),  i  =  0,  1,  — ,  N,  \j fa),  k  =  2,  3,  i  =  1,  2,  — ,  iV  },  (2.13) 

but  defer  more  specific  definitions  of  its  components  until  (2.16),  §3,  and  §4. 


Replacing  v  by  V  in  (2.2a)  and  subtracting  (2.12)  from  the  result  yields  the  ortho¬ 
gonality  relationship 

B(V,e)  :=  B(V,u-U )  =  0,  for  all  V  e  5^p.  (2.14) 

Replacing  v  by  e  in  (2.11)  and  subtracting  (2.14)  from  the  result  yields 

e{t)  -  B(G(t,-)  -  V,e),  for  all  V  e  Sq’p.  (2.15a) 

With  B(v,u)  continuous,  we  confine  t  to  %  and  obtain 

le(jcf)l  <  C\\G(Xi,-)  -  V'lli  \\e\\x,  for  all  V  e  5?’p,  i  =  0,  1,  -,  N.  (2.15b) 
The  estimate  (2.15b)  indicates  that  the  pointwise  error  of  the  Petrov-Galerkin  scheme 
(2.12)  can  be  reduced  by  making  \\e\\x  and/or  ||G (*,-,•)  -  V|li  small.  Satisfying  one  option 
usually  leads  to  the  other;  however,  trial  and  test  functions  can  differ  significantly  for 
singularly-perturbed  problems.  Using  (2.15),  Hemker  [13]  argued  that  Sq'p  should  be 
selected  to  produce  good  approximations  of  G  (*,•,*).  This  usually  implies  that  50,p 
should  accurately  represent  the  rapidly  varying  exponential  portion  of  G(xirx )  that  occurs 
at  x  =  xt  and  a  and/or  b  [16]. 


Example  2.2.  Consider  the  Petrov-Galerkin  solution  of  (2.1)  under  the  conditions  of 
Example  2.1.  Using  piecewise  linear  approximations  for  S% -1,  choose  a  basis  for  Sq’1 
that  exactly  satisfies  (2.10)  when  c  is  a  constant,  d  :=  0,  and  t  eUN,  i.e.,  choose 


Vt (£(*)),  if [*,•_!,*«) 

\| Vi(x)='  V_i(^)),  if  x  e  [Xi^i+i),  i  =  0,  1,  — ,  N , 
0,  otherwise 


(2.16a) 


where 


1  -  *-P(^)/2 

=  1  “  Vi(5).  Vi(5)  =  — ; - - — ■  (2.16b, c) 

1  -  e  w 

The  shape  functions  (2.16b,c)  reduce  to  the  linear  functions  (2.7c, d)  as  p  -4  0  and  have 
jump  discontinuities  at  £  =  -1  as  p  — »  °°  (cf.  Fig.  1). 

Using  (2.6,  7)  and  (2.16)  in  (2.12)  leads  to  the  discrete  system 

[1  +  -£-co(-£-)](cI+1  -  2c,  +  c,_j)  =  £-(c,+1  -  c,-^),  i  =  1,  2,  •••,  N  -  1,  (2.17a) 

with 

to(z )  =  cothz  -  —  (2.17b) 

Z 

and  the  boundary  conditions  (2.9a, b).  This  scheme,  which  is  identical  to  Il’in’s  [15]  finite 
difference  scheme,  yields  a  pointwise-exact  solution  of  (2.1)  under  the  conditions  of  this 
example. 

A  Petrov-Galerkin  scheme  such  as  (2.12)  would  require  a  major  re-coding  effort  to 
incorporate  into  a  state-of-the-art  adaptive  finite  element  software  system  (cf.,  e.g.,  Adjerid 
et  al.  [2]).  This  effort  would,  furthermore,  have  to  be  duplicated  for  different  singular  per¬ 
turbations.  Finite  element  schemes  typically  use  numerical  quadrature  to  evaluate  inner 
products;  thus,  instead  of  solving  (2.5),  a  solution  W*  e  S^  p  is  determined  from 

B*(V,W*)  =  (F,/)*,  for  all  V  e  S$-p ,  (2.18) 

where  the  *  indicates  that  integrals  are  evaluated  using  a  quadrature  rule. 

Example  2.3.  Hughes  [17]  recognized  that  the  diffusion  needed  to  stabilize  the 
piecewise-linear  Galerkin  solution  of  convection-diffusion  problems  could  be  added  by  a 
one-point  quadrature  rule  of  the  form 

l 

-1  L 


(2.19) 
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This,  when  used  with  (2.18),  yields  the  Petrov-Galerkin  discrete  system  (2.17).  The  qua¬ 
drature  rule  (2.19)  depends  on  p  and  approaches  the  midpoint  rule  (^  =  0)  in  the 
diffusion  limit  p  ->  0  and  the  Radau  formula  (^  =  sgnp)  in  the  convection  limit  P  -»  °°- 

Hughes’s  [17]  approach  only  requires  a  change  of  quadrature  rule  and  is,  thus, 
simpler  to  implement  than  a  Petrov-Galerkin  method.  Piecewise-polynomial  bases  (2.8a) 
are  used  for  both  trial  and  test  spaces  and  the  technique  may  be  extended  to  higher-order 
approximations  and  other  than  convection-diffusion  problems.  Thus,  we  seek  to  develop  a 
framework  for  the  design  of  quadrature  rules  that  are  appropnate  for  the  finite-element 
Galerkin  solution  of  singularly-perturbed  problems. 

3.  Quadrature-Based  Scheme  for  Convection-Diffusion  Problems.  Quadrature 
rules  for  convection-dominated  problems  are  developed  for  problems  having  the  form  of 
(2.1)  with  c(x)  *  0  and  d(x)  :=  0,  x  e  [, a,b ].  We  begin  by  defining  a  mapping  between 
the  polynomial  and  exponential  spaces  and  use  this  to  obtain  a  result  similar  to  (2.15) 
which  motivates  the  approach. 

DEFINITION  3.1.  Let  V(x)eS%’p  be  given  by  (2.6a),  with  c0  =  cN  =0  and  let 
F  :  Sq'p  -»  Sq’p  be  the  mapping 

V  -»  F(V)  =  V  =  qy/(x)  +  SI  cfytix).  (3.1) 

(=1  i=l k= 2 

LEMMA  3.1.  Let 

u‘  (x ) = £  c,V(* ) + £  i  c  *  >  <3-2a> 

i'=0  i=l  k=2 

be  an  element  of  Sg’p  that  satisfies  (2.18),  then 

e\x i)  :=  u(Xi)  -  U*(Xi )  =  B(G(Xi,-)  -  F(V),e*)  +  B*(V,U*)  -  B(F(V),U*)  + 

(F(V)J)  -  ( V,f  )*,  for  all  V  e  S^’p,  i  =  0,  1,  N.  (3.2b) 

Proof  Replacing  v  by  e*  in  (2.11),  adding  and  subtracting  B(F(V),e*)  to  the  right 
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side  of  the  result,  and  using  (2.18)  gives  (3.2b).  □ 

Using  (3.2b),  we  see  that  the  pointwise  error  can  be  reduced  by  (i)  selecting  'P  to  be 
a  good  approximation  of  the  Green’s  function  in  order  to  minimize  ||B(G(a-/,-)  -  V\\  and 
(ii)  designing  a  quadrature  rule  to  minimize  ||B*(V\W  )-B(F(V),U  )||  and 
\\(F(V),f)  -  (V,/)*||,  for  all  VeS%'p.  As  a  compromise,  we’ll  select  V  and  develop 
quadrature  rules  so  that  || B*(V ,W*)  —  B(F(V),U  )||  vanishes  for  locally  constant- 
coefficient  problems.  Thus,  using  (2.6c)  to  transform  integrals  on  (*,_!,*,)  to  (-1,1),  we 
require 

£*($*,£*)  =  B(\\fk£l),  k  =  -1,  1,  2,  •••,  p,  l  =0,  1,  ",p,  (3.3a) 

where 

l 

B(v,u)  =  J  [v'ffla'©  +  £v©u'©]<©  (3.3b) 

-1  2 

and  the  quadrature  rule  has  the  form 

1 

//©<© 

-l 

Elemental  indices  i  on  the  cell  Peclet  number  p  and  the  polynomial  degree  p  have  been 
omitted  for  convenience. 

Setting  p  =  1  and  using  (3.3)  with  \{/±1(£)  given  by  (2.16b, c)  and  $±1(£)  given  by 
(2.7b,c)  leads  to  Hughes’  [17]  one-point  quadrature  rule  (2.19).  Analysis  of  another  exam¬ 
ple  with  p  -  2  will  further  illustrate  the  technique  without  the  algebraic  complexity  of  the 
more  general  case  to  be  developed  subsequently. 

Example  3.1  ( Two-Point  Quadrature  Rule).  Bases  for  the  trial  and  test  spaces  for  a 
two-point  (n  =  2)  quadrature  rule  are  selected  as 

*={$-!.$!.  $2).  ^  V2  }• 


where 


(3.5a,b) 
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¥2(4)  = 


¥i(4)  ~  0i(4) 

“2 


(3.5c) 


and  the  scaling  factor  0C2  must  be  determined  (cf.  Fig.  1).  The  bases  (3.5a)  and  (3.5b) 
agree  except  for  the  last  component  which  contains  the  exponential  contribution  to  the 
Green’s  function.  The  shape  functions  <j>2  is  shown  in  Fi§-  L  As  p  0,  \j /2(4)  tends  to 
the  polynomial  shape  function  $2($)  and  as  p  -»  \jr2(§  becomes  linear  outside  of  an 
O  (1/p)  boundary  layer  at  ^  =  -1. 


Setting  p  =  2  in  (3.3)  and  using  (3.5)  yields  the  conditions 


W&i  +  =  J  4^4,  *=0,1,2, 

-1 


l  =  U2. 

This  nonlinear  system  is  readily  solved  to  obtain 

242  Ax 


(3.6a) 


(3.6b) 


3m  3  /  3m  3  _  ”  2q2 

«« -rW'T^i'  =  (3-7a‘d) 


where 


“2  =  Wlh*' 2X  m3  =  f( 


4'  1  -h 


(3.7e,f) 


®(p/2)  p 

and  co(z)  was  defined  by  (2.17b).  As  p  — >  0,  (3.7)  becomes  Gauss-Legendre  quadrature 
2  =  ±  1/V3,  Wj  =  W2  =  1)  and  as  p  -»  <*>,  (3.7)  becomes  Radau  quadrature 

(§!  =  -1/3,  ^2  =  1,  Wj  =  3/2,  W2  =  1/2). 

We  proceed  in  the  same  manner  with  a  general  ( n  -point)  quadrature  rule,  selecting 

O  =  {  0_1,  <j>i,  $2’  §n  }’  'P  =  { $-i>  $1*  $2>  — »  $n— 1>  ¥n  }*  (3.8a,b) 


where 


*  ...  ¥n-l(4)  -  $n-l(4) 

¥„(4)  = - -  n 


The  shape  functions  vjr3(£)  and  \i>4(^)  are  shown  in  Fig.  1. 


(3.8c) 
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The  use  of  (3.8)  with  (3.3,  4)  yields  the  system 


\  k  =0,  ",  2n  -  2, 

/=i  -l 

=  l  =  n  -  1,  «, 

which  may  be  written  in  the  more  explicit  form 


=  T 


1 4©4”-2<i5 

-1 


J  [Vi(S)-0i©]^>  if«  =  2 

-l 

J  V„-i(^”-2^  if»  >2 


(3.9a) 


(3.9b) 


(3.10a) 


£  W^2*-1  =  J 


(3.10b) 


/=1  C;(0)  -1 

The  nonlinear  system  (3.9a,  10)  may  be  solved  for  n  <4  using  a  computer  algebra 
system  such  as  MAPLE.  With  n  =  3,  for  example,  the  integration  points  l  =  1,  2,  3, 
are  determined  as  the  roots  of 


43_i^2_!^  +  i^  =  0, 


with 


3m  Wl5 


32  .  1  45  rj 

m5  =  - a3  -  - 

675  m3  2p  4 

When  n  =  4,  the  integration  points  are  the  roots  of 


105m  7  3 

+  £=0' 


with 


32  r  8  2  ^  „  - 

mn  =  t=t(- - a4  - 


15m5V35 


(3.11a) 


(3.11b,c) 


(3.12a) 


(3.12b, c) 
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While  the  quadrature  weights  Wt  and  evaluation  points  £/,  l  =  1,  2,  •••,  n,  cannot  be 
explicitly  determined  as  functions  of  p  for  n  >  4,  apparently  they  tend  to  the  Gauss- 


Legendre  weights  and  points  as  p  ->  0  and  to  the  Radau  weights  and  points  as  p  ->  ~. 
This  situation  is  favorable  to  adaptive  h-refinement  since  (as  we  shall  show  by  example) 
the  quadrature-based  finite  element  method  is  very  stable  in  the  convective  limit  and  accu- 

rate  in  the  diffusion  limit. 

3.1.  Computational  Results.  We  appraise  the  quadrature-based  finite  element 
method  (2.18)  by  applying  it  to  problems  involving  one  ordinary  and  two  partial 
differential  systems.  In  each  case,  the  Green’s  function  used  to  develop  the  quadrature 
rale  (3.4)  is  inaccurate;  hence,  we  hope  to  show  that  *  only  has  to  capture  the  essence  of 
the  singular  portion  of  the  adjoint  space  to  obtain  the  desired  stability.  We  also  show  that 
Radau  quadrature  provides  essentially  the  same  benefits  as  more  complex  quadrature  rules 
depending  on  the  cell  Peclet  number.  As  in  the  derivation,  the  local  degree  of  the  polyno¬ 
mial  approximation  p  is  the  same  as  the  number  of  points  n  used  for  the  quadrature  rule. 

In  each  example,  pointwise  errors  are  measured  in  the  discrete  maximum  norm 

le*  L  =  max  le*(*)l  (3-13) 

xenN{ 

where  Nx  is  the  number  of  elements  in  the  coarsest  mesh  used  to  solve  the  problem.  Spa¬ 
tial  complexity  is  measured  by  the  degrees  of  freedom  which,  with  Dirichlet  boundary 

data,  is  M  according  to  (2.8b). 

Example  3.2.  Consider  the  two-point  problem  [13] 

-e«"  -  xu  =  £TC2cos(7u:)  +  kx  sin(7U ),  -1  <  x  <  1,  (3.14a) 


w(-l)  =  —2,  k(1)  =  0, 


(3.14b,c) 


which  has  the  exact  solution 


u(x)  =  cos(tuc)  + 


erfpc  /VH) 
erf(l/V2e) ' 


(3.14d) 


The  solution  (3.14d)  is  smooth  outside  of  a  shock  layer  in  the  turning-point  region  near 


x  =  0. 
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We  solved  (3.14)  on  uniform  meshes  having  20,  40,  and  80  elements  with  piecewise 
polynomials  of  uniform  degrees  p  =  1,  2,  3,  4  and  values  of  e  =  10  2,  10-6,  and  10  . 

Maximum  errors  on  the  (Nx  =)  20-element  mesh  using  the  finite  element  method  with  the 
p-dependent  quadrature  rules  derived  from  (3.4,  8-10)  and  with  Radau  quadrature  rules  are 
presented  in  Table  1  as  functions  of  p,  N,  and  £.  In  Fig.  2,  we  display  \e  L  with 
e  =  10*6  as  a  function  of  the  degrees  of  freedom  for  p  ranging  from  1  to  4  using  (3.4,  8- 
10)  and  Radau  quadrature.  Similarly,  with  N  =  80,  we  display  \e  L  as  a  function  of  £ 
for  /?  =  1  to  4  in  Fig.  3.  Finally,  the  finite  element  solution  using  Radau  quadrature  is 
compared  with  the  exact  solution  when  £  =  ICT6 ,  N  =20,  and  p  =  1  to  4  in  Fig.  4. 


TABLE  1 

Maximum  errors  \e*L  for  Example  3.2  using  integration  formula  (3  4,  8-10) 
and  Radau  quadrature  on  uniform  N-element  meshes  with  piecewise  polynomials 

of  uniform  degree  p. 


Eas.  (3.4,  8-10) 

Radau 

- - - - - - =-= — hft — 

p 

N 

£  =  10"2 

£  =  10 

£  =  10^ 

E  =  10 

20 

0.109(  0) 

0.149(  0) 

0.1 82(  0) 

U.148(  U) 

1 

40 

0.402(  -1) 

0.765(  -1) 

0.98 1(  -1) 

0.764(  -1) 

80 

0.119(  -1) 

0.388(  -1) 

0.5 19(  -1) 

0.387(  -1) 

20 

0.1 24(  -2) 

0.142(  -3) 

0.526(  -2) 

0.142(  -3) 

40 

0.843(  -4) 

0.179(  -4) 

0.734(  -3) 

0.178(  -4) 

80 

0.533(  -5) 

0.224(  -5) 

0.988(  -4) 

0.223(  -5) 

20 

0.3 18(  -4) 

0.769(  -5) 

0.167(  -3) 

0.107(  -4) 

3 

40 

0.659(  -6) 

0.483(  -6) 

0.560(  -5) 

0.675(  -6) 

80 

0.1 18(  -7) 

0.302(  -7) 

0.188(  -6) 

0.422(  -7) 

20 

0.538(  -6) 

0.400(  -7) 

0.444(  -5) 

0.400(  -7) 

4 

40 

0.223(  -8) 

0.136(  -9) 

0.34 1(  -7) 

0.942(-ll) 

80 

0.875(-ll) 

0.302(-ll) 

0.272(  -9) 

0.186(-10) 

Finite  element  solutions  on  UNl  displayed  in  Table  1  and  Figs.  2  and  3  have  no 
spurious  oscillations  for  all  cell  Peclet  numbers.  Nodal  convergence  improves  as  p 
increases;  however,  there  is  an  0(e)  error  that  cannot  be  removed  without  proper  resolu¬ 
tion  of  the  solution  within  the  turning-point  region.  (This  phenomena  also  occurs  with 
boundary  layer  problems.)  The  special  quadrature  rules  (3.4,  8-10)  produce  solutions  that 
are  slightly  better  than  those  obtained  with  Radau  quadrature,  but  the  difference  may  not 
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FIG.  2.  Maximum  errors  \e*\„  vs.  degrees  of  freedom  M  for  Example  3.2  using  in¬ 
tegration  formula  (3.4,  8-10)  (left)  and  Radau  quadrature  (right).  Uniform  mesh  and  order 
computations  correspond  to  p  =  1  (O),  2  (A),  3  (+),  and  4  (x). 

be  worth  the  added  expense.  Results  presented  in  Fig.  4  show  that  the  finite  element- 
Radau  solution  has  some  excess  diffusion  when  p  =  1  and  some  spurious  oscillations 
when  p  >  1;  however,  these  undesirable  effects  are  confined  to  the  two  elements  contain¬ 
ing  the  turning  point.  The  oscillations  decrease  in  magnitude  as  p  increases  and  the  poly¬ 
nomial  basis  provides  a  better  approximation  to  the  exponential  boundary  layer  behavior. 
Approximations  would  likewise  improve  were  the  mesh  refined  in  the  turning-point  region. 
Global  accuracy  away  from  the  turning  point  is  very  high. 

Example  3.3.  In  order  to  appraise  the  method’s  suitability  for  use  with  transient 
problems,  consider  Burgers’  equation 


t 


1/e  1/e 


FlG.  3.  Maximum  errors  \e*\oa  vs.  l/£  for  Example  3.2  using  integration  formula 
(3.4,  8-10)  (left)  and  Radau  quadrature  (right).  Uniform  mesh  computations  with  N  =  80 
were  performed  with  uniform  degrees  of  p  =  1  (O),  2  (A),  3  (+)  and  4  (x). 

ut+uux=euxx,  0  <  x  <  1,  t  >  0,  (3.15a) 

with  initial  and  Dirichlet  boundary  conditions  prescribed  so  that  the  exact  solution  is  the 
traveling  wave  [18] 


u(x,t)  = 


0.le~A  +  0.5e~B  +  e~c 


e~A  +  e  B  +  e  c 


(3.15b) 


with  e  =  3  x  10  3  and 


A  =M1(X_  o.5 +4.950,  B=—  (x- 0.5 +  0.750,  C  =  — (x  -  0.375) 


0.25 


0.5 


C 
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FIG  4  Finite  element  solution  using  Radau  quadrature  and  exact  solution  of  Exam¬ 
ple  3.2  with  e  =  KT6,  N  =  20  and  p  =  1  (O),  2  (*),  3  (+)  and  4  (x). 

The  cell  Peclet  number,  needed  for  use  with  (3.4,  8-10),  was  determined  from  (2.9d) 
with  c  replaced  by  the  average  of  the  convective  velocity  U*  at  both  ends  of  an  element. 
We  solved  (3.15)  for  0  <  t  <  0.5  on  a  uniform  80-element  mesh  using  piecewise  polyno¬ 
mials  having  uniform  degrees  of  one  to  four.  Temporal  integration  used  the  backward 
difference  code  DASSL  [19]  with  an  error  tolerance  of  10-10.  The  exact  solution  and  the 


r 
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finite  element  solution  using  (3.4,  8-10)  are  compared  as  functions  of  x  at  r  =  0.5  in  Fig. 
5  for  p  =  1,  2.  Results  with  p  =  3  and  4  could  not  be  distinguished  from  the  exact  solu¬ 
tion.  The  piecewise  linear  solution  shown  on  the  left  of  Fig.  5  has  too  much  dissipation 
and  an  incorrect  wave  speed.  The  piecewise  quadratic  solution  has  remedied  these 

deficiencies. 


FIG  5  Finite  element  (dashed)  and  exact  (solid)  solutions  of  Example  3.3  at  t  =  0  5 
using  a  uniform  80-element  mesh  and  piecewise  polynomials  of  uniform  degree  one  (left) 

and  two  (right). 


Example  3.4.  Consider  the  nonlinear  initial-boundary  value  problem  [20] 


£Uf  +  u(u2-l)ux  +  u  = 

0  <  x  <  1,  t  >  0, 

(3.16a) 

u  (x  ,0)  =  (x  +  3)/2, 

0  <  x  <  1, 

(3.16b) 

w  (0,0  =  3/2,  u  (1,0 

=  2,  t  >  0. 

(3.16c,d) 

The  solution  of  (3.16)  tends  to  a  steady  state  as  t  -*  ~  and  time  has  been  scaled  to 
quickly  reach  this  limit.  An  analysis  when  0  <  e  <sc  1  [20]  reveals  that  the  steady  solution 
features  an  interior  layer  near  x  =  0.096,  a  comer  layer  near  x  =  0.333,  and  a  boundary 


r 
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layer  at  .v  =  1. 

Equations  (3.16)  were  solved  for  t  e  (0.0.1]  with  £  =  10-3  using  adaptive  h- 
refinement  with  uniform  piecewise  quadratic  polynomials  [2],  an  Hx  spatial  error  tolerance 
of  0.1,  and  a  temporal  error  tolerance  of  10-4.  The  finite  element  solution,  computed  with 
(3.4.  8-10)  using  cell  Peclet  numbers  as  specified  in  Example  3.3,  and  the  mesh  con¬ 
structed  by  the  adaptive  procedure  are  shown  in  Fig.  6.  The  solution  differed  from  one 
computed  using  Gauss-Legendre  quadrature  by  less  than  2x10  3  in  strain  energy,  how¬ 
ever,  the  solution  obtained  using  (3.4,  8-10)  used  20%  fewer  space-time  degrees  of  free¬ 
dom.  The  adaptive  mesh  is  concentrated  in  layers  and  little  element  removal  was  neces¬ 
sary.  Meshes  track  evolving  layers  from  the  smooth  initial  data.  The  solution  obtained 
with  Gauss-Legendre  quadrature  performed  well  here  because  h-refinement  occurred 
quickly  as  layers  developed. 


FIG.  6.  Finite  element  solution  of  Example  3.4  at  t  —  0.2  using  the  quadrature  rules 
(3.4,  8-10)  with  adaptive  h-refinement  and  p  =  2  (left)  and  the  mesh  used  for  the  adaptive 
computation  (right). 
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4.  Quadrature-Based  Scheme  for  Reaction-Diffusion  Problems.  We  develop  qua¬ 
drature  rules  for  reaction-diffusion  problems  of  the  form  of  (2.1)  with  c(x)  :=  0  and 
d(x)  >  0  using  the  approach  described  in  §3.  The  analog  of  Lemma  3.1  has  a  slightly 

more  complex  form. 


LEMMA  4.1.  Let  d>  of  (2.8a)  and 

H  :=  { ti/Cx),  i  =0,  1,  •••,  N,  r\t(x),  k  =  2,  3,  J  =  L  2,  —,N  },  (4.1) 


(F(V)J)  -  for  all  V  e  Sq,p  ,  i  =0,1,  •••,#.  (4.3) 

Remark.  The  elements  of  H  will  be  chosen  to  be  identical  to  those  of  <b  except  for 

nf-  which  will  differ  from  <i>f\  pt  >  1,  by  a  scaling  factor,  i  =  1,  2,  N.  The  scaling  is 
necessary  to  avoid  a  constant  difference  between  exact  and  numerical  inner  products. 

Proof.  A  direct  computation  following  the  steps  of  Lemma  1  yields  the  result.  □ 

As  a  function  of  jc,  the  Green’s  function  in  this  case  has  0(<e)  boundary  layers  on 
both  sides  of  jc  -t  and  becomes  unbounded  as  0(1/Ve)  as  e  ->  0  [16].  Requiring 
\\B*(V  W*)  -  B(F(V),U*)\\  to  vanish  for  locally  constant-coefficient  problems  and 

transforming  via  (2.7a)  to  the  canonical  element  (—1,1),  we  obtain 

S*(0fe.^)  =  5(¥fc.fl/)>  k'1  =-1’  !»2»  "  ’P' 


(4.4a) 


where  now 


B (v  ,u )  =  J  [v'  (Qu  (£)  +  (£)m  (^)]  d%,  a  =  h *\pt.  (4.4b, c) 

-l  4 

For  this  self-adjoint  problem,  we  select  a  symmetric  quadrature  rule  of  the  form 


J  f<S)d% 


=  Wq/^q)  + 


int  (n  12) 

I  wk[f<^) +  /<&)] 


k=  1 


(4.5) 


where  0  =  <  ^n,(n/2)  ^  1  md  int^  denotes  the  inte§er  VMof  x-  The  mle 

is  to  have  precisely  n  points,  so  Wq  is  zero  when  n  is  even. 


Being  unable  to  obtain  results  of  the  generality  of  those  in  §3,  we  examine  specific 
cases  with  p  =  1,  3,  and  5.  Quadrature  rules  having  the  form  of  (4.5)  were  either  incom¬ 
patible  with  (4.4)  or  produced  integration  points  outside  of  [-1,1]  with  even  values  of  p. 

4.1.  Piecewise-Linear  Approximation.  With  a  linear  polynomial,  the  bases  for  the 
trial  and  test  spaces  are  selected  as 

4>  =  H  =  {$_!,  },  ❖=  (V_i,  ViJ,  (4.6a, b) 

where  (j)*.,  k  =  -1,  1,  are  given  by  (2.7b, c)  and 

-  ,  =  sinho(  1  ±  ^)/2  (4.6c) 


Substituting  (4.6)  into  (4.4a)  and  using  (4.5)  yields  the  independent  relations 

Xo  =  —  tanh(c/2),  X2  =  —  [coth(c/2)  -  -^-tanh(a/2)] 

G  O  Gl 


(4.7) 


where 


X,  = 


-1 


int  (n  12) 

=  wVjo+  £  wiK-V  +&)']■ 

*=1 


(4.8) 


Combined  use  of  (4.8)  and  (4.7)  indicates  that  there  is  no  quadrature  rule  with  n  =  1  satis¬ 
fying  (4.7).  A  two-point  quadrature  rule  exists  with  W0  =  0  and 
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W  j  =  1,  ^  =  Vcoth2a/2  -  4/a2.  (4.9) 

(Scaling  of  is  arbitrary  and  we  choose  unity  for  convenience.)  In  the  limits  as  a  tends 
to  zero  and  infinity,  ^  approaches  V2/3  and  1,  respectively.  The  former  case  produces  an 
order  one  quadrature  rule  and  the  latter  case  yields  the  trapezoidal  rule  (or  the  two-point 
Lobatto  formula). 


4.2.  Piecewise-Cubic  Approximation.  The  bases  of  the  approximation  spaces  for 
cubic  trial  functions  are  selected  as 


6  :=  {$_1,  $1.  $2,  fe).  H:=  {$-1.  $1.  $2’  a3<t>3  }.  ^  :=  {$-!’$!’  V2>  V3}>(4-10) 


where  the  scaling  a  3  is  to  be  determined.  Continuity  considerations  require  \j/,  (±l)  =  0, 
i  =  2,  3.  Additionally,  <j>f(£),  i  >  2,  is  an  odd  or  even  function  of  £,  when  i  is,  respec¬ 
tively,  odd  or  even  (cf.  (2.7d)).  These  restrictions  and  the  use  of  (4.4a)  with  k  =  2, 
/  =  -1,  1,  3,  and  k  =  3,  /  =  -1,  1,  2,  imply  that  y2  is  even  ¥3 is  odd-  Thus’ 


P3, 


Vz(^)  =  Pzt1  "  ¥i(S)  -  V-i(0].  V3(S)  =  PsK  -  ¥i®  +  ¥-i(4)]- 
Imposing  conditions  (4.4a)  and  solving  the  resulting  system  for  X0,  X2,  P2,  a3> 
and  X6  while  using  (4.11)  and 


“¥"±1  +  "^"¥±1  -  (4.12) 

we  find 

1  8 

X0  =  2,  X2  =  2/3,  P2  = - T - >  a3  =  5(1  -  P2  +  — ), 

V6(l  -  j  \jjr1(^) z/ 4) 

-l 


(4.13a-d) 


X4  =  m4  =  -|(2p2-  1  -  ^),  p3--i  1 


a3 


(4.13e,f) 


1  -  3  J  ^\j>i (£)<*£ 

-l 


X6  =  m6  =  -|[-1  +  +  3(1  -  -^)m4  -  4yjjM. 


(4.13g) 


(4. 14a,b) 


We  may  use  (4.8)  to  write  (4.13a,b,e,g)  in  the  more  explicit  form 


W, 


int(n/2) 

+  2  %  Wk 

k=i 


=  2, 


in!  (n  12) 

£ 

*=i 


w. 


int  (n  12) 

2  = 


k=\ 


m4 

T’ 


inr  (n 12)  m  e. 

I  wkti  =  ~Y- 

k= i  z 


(4.14c,d) 


A  quadrature  rule  satisfying  (4.14)  with  n  =  4  and  W0  =  0  exists  but  has  ^2  outside  [-1,1] 
and,  hence,  is  useless  for  our  application.  When  n  =  5,  (4.14)  describe  a  one-parameter 
family  of  quadrature  rules  that  we  specify  by  selecting  £2  =  L  In  this  case’ 


„  „  ,2  m4  ~  m6 

to  -  °-  tl  -  2/3  _  m4  .  $2  • 


(4.15) 


W0  =  2(1  -  Wt  -  W2),  Wj  = 


1/3  -  W2 

— tT" 


,  W2  = 


3m 4  -  2 m6 
6m  4  —  3m  6  +  2 


With  all  integration  points  within  [-1,1],  we  find  the  order-five  formula 


$o  =  0,  4,  =  ^=.  5j=1.  w0  =  -8.  w,= 


51,  2  51 


(4.16) 


(4.17) 


as 


a  -4  0.  As  cr  — >  °°,  (4.15,  16)  give 


5o  =  0,  =  fe-1.  W0  =  0,  W,  =  1,  Wj-O  (4.18) 

which  is  the  reduced  two-point  Gauss-Legendre  formula.  A  seven-point  quadrature  rule 
for  quintic  polynomials  is  not  reproduced  here  due  to  its  algebraic  complexity  [21].  It  has 
properties  that  are  similar  to  (4.15,  16).  For  example,  it  approaches  a  reduced  four-point 
Gauss-Legendre  formula  as  a  -»  00 . 

These  surprising  results  are  difficult  to  understand.  Based  on  the  analysis  of  §3  and 
§4.1,  we  would  have  expected  to  find  formulas  that  approach  Gauss-Legendre  quadrature 
as  C  — >  0  and  Lobatto  quadrature  as  a  — >  Indeed,  we  shall  show  in  §4.3  that  Lobatto 

quadrature  rules  produce  stable  and  accurate  results  for  all  values  of  <7.  However,  unlike 
the  Radau  rules  discussed  in  §3,  the  Lobatto  rules  apparently  do  not  follow  from  the 


-  25  - 


formalism  of  §4.1  and  here. 

4.3.  Computational  Results.  Finite  element  solutions  are  compared  using  the  special 
formulas  derived  in  §4.1  or  §4.2  and  Lobatto  quadrature.  The  Lobatto  formulas  are 
chosen  to  have  p  +  1  points  for  piecewise  polynomial  approximations  of  degree  p . 


Example  4.1.  Consider  the  two-point  problem  [16] 


-eV'  +  (x2  +  £)«  =  (. x 2  +  e)(l  +  sinnx)  +  eVsimtx,  0  <  x  <  1, 

(4.19a) 

n(0)  =  m(1)  =  0, 

(4.19b,c) 

which  has  the  exact  solution 

u (pc')  =  1  +  sinrcc  -  uH(x) 

(4.19d) 

where 

1  {(1  e-'^W(xl<i)e^l2‘  +  n-e-'l7zW(<m]e<'- 

HK  '  erf(Vl/e) 

-x2)I2e  j 

(4.19e) 

and 

W(z)  =  ez2erfcz. 

(4.19f) 

This  solution  features  boundary  layers  at  x  =  0  and  1.  The  boundary  at  x  -  0  is  nearly  a 


second-order  turning  point. 

We  solved  (4.19)  using  the  finite  element  method  with  the  special  quadrature  rules 
(4.5)  on  uniform  meshes  having  10,  20,  40,  and  80  elements  with  piecewise  polynomials 
of  uniform  degrees  p  =  1,3,  5.  Maximum  errors  on  the  10-element  mesh  are  presented 
for  e  =  10'3,  10-5,  and  10-7  in  Table  2.  Similar  results  appear  in  Table  3  for  computa¬ 
tions  performed  with  Lobatto  quadrature.  A  display  of  le*L  as  a  function  of  the  degrees 
of  freedom  and  e  is  presented  in  Fig.  7  for  the  finite  element-Lobatto  solution.  Results 
with  the  special  quadrature  rule  (4.5)  had  a  similar  behavior.  Finally,  the  exact  and  finite 
element-Lobatto  solutions  with  e  =  10-5,  N  =  10,  and  p  =  1  to  4  are  compared  in  Fig.  8. 


Table  2 


Maximum  errors  \e*\oa  for  Example  4.1  using  integration  formula  (4.5)  on  uni¬ 
form  N-element  meshes  with  piecewise  polynomials  of  uniform  degree  p. 


p 

N 

1 

o 

II 

CO 

10"5 

10"7 

10 

0.269(  -2 r 

0.486(  -6)' 

0.487(-i(J) 

20 

0.37 1(  -3) 

0.143(  -7) 

0.143(-1 1) 

1 

40 

0.254(  -3) 

0.149(  -7) 

0.150(-1 1) 

80 

0.85 1(  -4) 

0.151(  -7) 

0.151(-11) 

10 

0.437(  -2) 

0.477(  -3) 

0.488(  -5) 

20 

0.382(  -4) 

0.339(  -6) 

0.1 37(  -9) 

3 

40 

0.272(  -6) 

0.1 11(  -8) 

0.519(-10) 

80 

0.127(  -7) 

0.1 32(  -9) 

0.3 19(- 10) 

10 

0.386(  -4) 

0.905(  -3) 

0.976(  -5) 

C 

20 

0.172(  -6) 

0.926(  -6) 

0.1 50(  -9) 

D 

40 

0.172(-10) 

0.401(-12) 

0.529(-10) 

80 

0.972(-13) 

0.597(- 12) 

0. 195(- 10) 

Table  3 

Maximum  errors  \e*  L  for  Example  4.1  using  Lobatto  quadrature  on  uniform 
N-element  meshes  with  piecewise  polynomials  of  uniform  degree  p. 


P 

N 

e  =  10~3 

10~3 

10~7 

10 

0.777C-2) 

0.998(  -6) 

0.999R0) 

1 

20 

0.201(  -2) 

0.297(-10) 

0.621  (-14) 

1 

40 

0.546(  -3) 

0.156(-10) 

0.155(-14) 

80 

0.140(  -3) 

0.391(-11) 

0.666(-15) 

10 

0.322(  -2) 

0.998(  -6) 

0.999(-10) 

20 

0.475 (  -3) 

0.947(-10) 

0.3 10(- 14) 

z 

40 

0.468(  -4) 

0.782(- 11) 

0.3 10(- 14) 

80 

0.361(  -5) 

0.195(-1 1) 

0.222(-14) 

10 

0.121(  -3) 

0.998(  -6) 

0.100(  -9) 

o 

20 

0.1 09(  -4) 

0.628(-10) 

0.488(-14) 

5 

40 

0.38 1(  -7) 

0.444(-14) 

0.621  (-14) 

80 

0.702(  -9) 

0.577(- 14) 

0.288(- 14) 

10 

0.1 18(  -3) 

0.995(  -6) 

0.1 00(  -9) 

A 

20 

0.636(  -6) 

0.602(- 10) 

0.999(-14) 

4 

40 

0.769(  -8) 

0.888(- 14) 

0.643(-14) 

80 

0.328(-10) 

0.1 13(-13) 

0. 1 15(- 13) 

10 

0.1 13(  -3) 

0.985(  -6) 

0.999(-10) 

r 

20 

0.1 13(  -6) 

0.528(- 10) 

0.106(-13) 

J 

40 

0. 156(- 10) 

0.999(-14) 

0. 177(- 1 3) 

80 

0.142(-13) 

0.126(-13) 

0.577(-14) 

-  27  - 


FIG  7  Maximum  errors  le*L  vs-  degrees  of  freedom  M  (left)  and  vs.  1/e  (right) 
for  Example  4.1  using  Lobatto  quadrature.  Uniform  mesh  and  order  computations 
correspond  to  p  —  1  (O),  2  (A),  3  (+),  4  (x),  and  5  (CH). 

As  with  convection-diffusion  systems,  solutions  on  nVl,  presented  in  Tables  2  and  3 
and  Fig.  7,  have  no  spurious  oscillations  for  all  values  of  c.  Results  using  Lobatto  qua¬ 
drature  are  either  comparable  or  superior  to  those  obtained  by  the  special  quadrature  rules 
derived  from  (4.5);  thus,  there  appears  to  be  little  advantage  of  using  the  more  complex 
a-dependent  rules.  Results  presented  in  Fig.  8  show  that  boundary  layer  errors  are 
confined  to  one  element  when  a  is  large.  The  piecewise-linear  solution  has  no  oscillations 
but  higher-order  solutions  have  spurious  oscillations  within  layers  that  decay  in  amplitude 


with  increasing  p . 


p-1  P=2 


p  =  3  P  =  4 


FIG.  8.  Finite  element  solution  using  Lobatto  quadrature  and  exact  solution  of  Exam¬ 
ple  4.1  with  e  =  10-5,  N  =  10  and  p  =  1  (O),  2  (*),  3  (+)  and  4  (x). 

5.  Two-Dimensional  Problems.  To  have  maximal  impact,  the  specialized,  Radau, 
and  Lobatto  quadrature  rules  developed  and  described  in  §3  and  §4  should  be  applicable 
to  multi-dimensional  transient  and  steady  singularly-perturbed  partial  differential  systems. 
We  appraise  their  suitability  in  this  regard  by  applying  a  tensor  product  of  the  one¬ 
dimensional  quadrature  rules  to  three  two-dimensional  elliptic  problems. 


r 
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Example  5.1.  Consider  the  convection-diffusion  equation 

-eAw  +  2ux  +  uy  =  /(x,  y),  (x,  y)e  (-1,1) x (-1,1),  (5.1a) 

with  f(x,y)  and  the  Dirichlet  boundary  conditions  selected  so  that  the  exact  solution  is 

u(x,y)  =  (1  -  g-^-^Xl  -  e-(1_y)/e)cos7t(x  +  y).  (5 -lb) 

This  solution  features  0(e)  boundary  layers  at  x  =  1  and  y  =  1. 

We  solved  (5.1)  with  e  =  10~3  and  1CT6  using  the  finite  element  method  with  a  tensor 
product  of  the  quadrature  rules  (3.4,  8-10)  on  uniform  square  meshes  having  8,  16,  and  32 
square  elements  per  edge  and  piecewise  bi-polynomial  approximations  of  uniform  degrees 
one  to  four.  The  maximum  pointwise  errors  measured  on  the  coarse  mesh  are  presented 
in  Fig.  9.  Solutions  are  computed  without  oscillations  for  all  combinations  of  e  and  N. 

Example  5.2.  Again  consider  (5.1a)  with/(x,y)  =  0  and  the  boundary  conditions 

u(x,0)  =  1,  w(x,l)  =  2,  0  <  x  <  1,  u  (0,y )  =  2,  u(hy)  =  1,  0<y<l.  (5.2) 
When  e  is  small  relative  to  unity,  the  solution  features  a  sharp  wave  front  that  propagates 
across  the  domain  at  an  angle  of  approximately  27°  with  respect  to  the  positive  x-axis. 

We  solved  (5.1a,  2)  with  £  =  10"3  using  the  tensor-product  quadrature  rules  (3.4,  8- 
10)  on  a  20x20  uniform  mesh  and  piecewise  bi-polynomial  approximations  having 
degrees  one  through  four.  Solutions  are  displayed  in  Fig.  10.  Like  Brooks  and  Hughes 
[22],  we  find  solutions  with  piecewise  bilinear  approximations  to  be  overly  diffusive. 
Higher-order  solutions  have  less  diffusion,  but  have  some  spurious  oscillations  near  the 
wavefront  that  decrease  in  amplitude  with  increasing  p.  Streamwise  upwinding  [22]  has 
been  used  with  low-order  approximations  to  remove  excessive  diffusion  near  fronts. 
Perhaps  a  similar  procedure  could  be  developed  to  further  reduce  the  oscillations  associ¬ 
ated  with  higher-order  approximations. 

Example  5.3.  Holland  [23]  suggested  the  resonance  problem 

0,  -2  <  x  <  1,  -3  <  y  <  3. 


-eA u  +  xux  +  yuy  = 


(5.3a) 


M 


M 


RG.  9.  Maximum  errors  l<?*L  vs.  degrees  of  freedom  M  for  Example  5.1  with 
E  =  10-3  (left)  and  10"6  (right)  using  a  tensor  product  of  the  quadrature  rules  (3.4,  8-10) 
on  uniform  square  meshes  with  piecewise  polynomials  of  uniform  degree  p  =  1  (O),  2 
(A),  3  (+),  and  4  (x). 

u  (x  ,-3)  =  3,  u(x, 3)  =  5,  -2  <  x  <  1,  (5.3b,c) 

u(l,y)  =  4,  u(-2,y)  =  6,  -3  <  y  <  3.  (5.3d, e) 

The  usual  singular-perturbation  theory  would  indicate  that  the  solution  of  (5.3)  has  boun¬ 
dary  layers  near  the  edges  and  is  constant  in  the  interior  of  the  domain;  however,  this 
theory  cannot  determine  the  constant’s  value.  Grassman  and  Matkowsky  [24]  used  a  vari¬ 
ational  approach  to  determine  the  unknown  constant  as  a  weighted  average  of  the  boun¬ 
dary  data  at  points  that  are  closest  to  the  origin.  For  the  prescribed  data  and  small  values 
of  e,  there  will  be  boundary  layers  except  near  (1,0),  which  is  the  closest  point  to  the  ori- 


o 


gin.  The  solution  in  the  interior  of  the  domain  will  asymptotically  be  u(l,0)  =  4. 

Holland  [23]  suggested  this  problem  to  test  numerical  techniques  because  of  its 
numerous  computational  difficulties.  A  transient  embedding  of  (5.3a)  in  a  parabolic  prob¬ 
lem  [2]  converges  at  an  exponentially  slow  rate  in  e.  When  solving  the  steady  problem 
(5.3),  this  difficulty  gives  rise  to  ill-conditioned  discrete  systems  for  large  cell-Peclet 
numbers.  Thus,  direct  solution  techniques  will  be  sensitive  to  round-off  error  accumula¬ 
tion  and  iterative  solution  strategies  will  quickly  converge  to  a  constant  interior  solution, 
but  take  exponentially  long  to  find  the  correct  value. 

We  solved  (5.3)  with  e  =  10"3,  using  adaptive  h-refinement  [2]  with  piecewise  bi¬ 
polynomials  of  degrees  one  through  four.  The  solution  with  p  =  1  is  shown  in  Fig.  11. 
Solutions  did  not  display  any  spurious  oscillations. 

6.  Discussion.  We  have  developed  a  framework  for  applying  the  finite  element 
method  with  high-order  approximations  to  singularly-perturbed  differential  systems.  The 
method  utilizes  symbolic  techniques  to  construct  quadrature  rules  for  a  class  of  singular- 
perturbation  problems  and,  herein,  we  consider  convection-diffusion  and  reaction-diffusion 
systems.  Quadrature  rules  for  convection-diffusion  systems  tend,  as  expected,  to  Gauss- 
Legendre  and  Radau  integration,  respectively,  as  the  cell-Peclet  number  tends  to  zero  and 
infinity.  Quadrature  rules  are  less  understood  for  reaction-diffusion  systems.  Formulas  for 
odd-degree  polynomial  approximations  produced  stable  results  but  appeared  to  be  sub- 
optimal.  Formulas  for  even-degree  polynomials  had  evaluation  points  off  the  element. 

Use  of  Radau  and  Lobatto  quadrature  rules  worked  extremely  well  for,  respectively, 
convection-  and  reaction-diffusion  problems.  Large  errors  were  confined  to  elements  con¬ 
taining  boundary  or  interior  layers  for  all  meshes,  orders,  and  singular-perturbation  param¬ 
eters  tested.  This  is  in  contrast  to  the  Radau-  and  Lobatto-based  collocation  methods  of 
Ascher  and  Weiss  [25]  who  found  oscillations  when  boundary  layers  were  not  adequately 
resolved.  Furthermore,  nodal  convergence  of  the  Radau-  and  Lobatto-based  finite  element 
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FIG.  11.  Solution  of  Example  5.3  with  e  =  10  3  using  adaptive  h-refinement  on  an 
initial  10  x  10  uniform  mesh  with  piecewise  bilinear  approximations. 

procedures  appears  to  be,  respectively,  at  rates  of  h2p~l  and  h2p  in  the  diffusion  limit. 
Thus,  nodal  superconvergence  would  seem  to  be  present  for  both  quadrature  rules  when 
p  >2.  Global  rates  of  hp  are  optimal  in  energy  in  the  diffusion  limit  but  are,  as  yet,  unk¬ 
nown  in  the  singularly-perturbed  limit.  Observed  accuracy  is  so  high  in  this  case  that  esti¬ 
mation  is  not  possible.  There  is  little  apparent  advantage  to  using  the  more  complex 
integration  procedures  of  §3  and  §4,  especially  with  reaction-diffusion  problems.  The  qua¬ 
drature  framework  (§3  and  §4)  could,  of  course,  be  useful  with  other  singularly-perturbed 
problems.  It  additionally  provides  insight  as  to  why  Radau  quadrature  is  successful  with 
convection-diffusion  problems. 


Numerical  evidence  suggests  that  the  quadrature-based  methods  work  more  robustly 
than  anticipated.  Thus,  for  example,  methods  based  on  a  singular-perturbation  analysis  of 
a  constant-coefficient  two-point  boundary  value  problem  provide  stable  and  accurate  solu¬ 
tions  of  two-point  problems  involving  turning  points  and  of  partial  differential  equations. 
Having  stable  high-order  methods,  it  becomes  possible  to  use  efficient  adaptive  hp- 
refinement  procedures  and  we  intend  to  investigate  this  possibility.  Several  aspects  of  the 
approach  are  in  need  of  additional  analysis  before  this  can  be  done.  A  posteriori  error 
estimates,  used  to  guide  adaptive  enrichment,  are  needed  for  each  method  and  quadrature 
rule.  It  is  likely  that  such  estimates  can  be  developed  by  p-refinement  procedures  [10]. 
Indeed,  Biswas  et  al.  [26]  used  a  Radau  polynomial  to  construct  error  estimation  formulas 
for  hyperbolic  conservation  laws. 

A  priori  error  estimates  are  also  needed  for  the  various  methods  in  the  different 
parameter  regimes  characterized  by,  e.g.,  the  cell-Peclet  number.  Streamwise  upwinding 
should  also  be  investigated  as  a  possibility  for  reducing  oscillations  near  nonuniformities 
that  are  oblique  to  the  computational  mesh  (cf.  Example  5.2).  Developing  such  formulas 
for  high-order  approximations  could  be  a  challenging  proposition  since  streamline  curva¬ 
ture  may  be  necessary  to  maintain  high  order. 
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